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ABSTRACT 

We use a multi-dimensional hydrodynamics code to study the gravitational inter- 
action between an embedded planet and a protoplanetary disk with emphasis on the 
generation of vortensity (potential vorticity) through a Baroclinic Instability and sub- 
sequent development of Rossby- Wave-Instability (RWI). It is found that the generation 
of potential vorticity is very common and effective in non-barotropic disks through the 
Baroclinic Instability, es pecially within the coorbital region. Our results also comple- 



ment previous studies by lKoller et alJ (120031 ) that non-axisymmetric RWIs are likely to 



develop at local minima of potential vorticity distribution that are generated by the 
interaction between a planet and a inviscid barotropic disk. This second instability ap- 
pears to be very common and robust, regardless of the equation of state, initial density 
distribution, and rotational law of the disk. The development of RWIs results in non- 
axisymmetric density blobs, which exert stronger torques onto the planet when they 
travel in the vicinity of the planet. As a result of that, large amplitude oscillations are 
introduced to the time behavior of the total torque acted on the planet by the disk. 
In our current simulations, RWIs do not change the overall picture of inward orbital 
migration but bring in a non-monotonic behavior to the migration speed. As a side 
effect, RWIs also introduce interesting structures into the disk. These structures may 
help the formation of Earth-like planets in the Habitable Zone or Hot Earths interior 
to a close-in giant planet. 

Subject headings: accretion, accretion disks — Baroclinic Instability — Rossby-wave in- 
stability — hydrodynamics — numerical methods — planetary systems: protoplanetary 
disks 
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Introduction 



Although theorists did not reahze its relevance to the origin of planet when the study on the 
gravit ational interaction between a protoplanetary disk and an embedded planet just began in late 
1970s (jGoldreich Tremainelll978l ). it is now one of the key ingredi ents of current standard theor y 
of planet formation. According to standard core-accretion theory (|Safronovlll969l : lLissaueiill993l ). 
the formation of planets in circumstellar disks around T Tauri stars consists of the formation of 
planetesimals via collisions/coalitions of dust grains in the early stage and then gravitationally 
accretion after they accumulate enough mass. The life time of this T Tauri star phase is estimated 
to be short (< 10^ years). In order for cores of protoplanets with masses comparable to that of 
Jupiter to accumulate enough mass in the T Tauri stage, it is thought that their cores have to form 
outside the so-called "ice lin e", where the distance to the central star is large (typically beyond 
~ 4 AU, see llda &: LinI |2004| ) so that the temperature is low enough to allow the condensation of 
gas materials to solid ice. These additional solid dusts helps to increase the dust coagulation speed 
and shorten the time needed to form cores of protoplanets. 



Since the discovery of the first Jupiter-mass planet (iMayor &: Queloall995l ) orbiting the solar- 
type star 51 P eg, more than 200 extrasolar planets have been discovered around the nearby stars 
within 200 pc (jButler et al.ll2006l . The Extrasolar Planets Encyclopaedia). The diversity exhibited 
by these planetary systems shows that the masses of these planets range from Jupiter-mass to 
Neptune-mass. The most prominent characteristics is that a large number of the host stars are 
surrounded by the so-called "hot- Jupiters" and close-in super-Earths. Around 80% of the extrasolar 
planets are in orbits with semi-major axes in the range 0.01 ^ a < 2.5 AU, and ~ 25% of the total 
population are short-period planet^ with a < 0.1 AU. This has brought one of the most interesting 
puzzles to theorists: if protoplanets had formed in a disk region beyond ~ 4 AU from the central 
star, how did the observed extrasolar planets end up with orbits that are so close to their host 
stars? If the standard theory for the formation of protoplanetary cores holds (Pollack et al. 1996; 
Ida & Lin 2004), then the giant pl anets (like hot- J upiters) must have undergone inward orbital 
migration to their current locations (jLin et al.lll996l ). 



S everal mechaiiisms have been propose d as the drive for the migration. Some authors (jWeidenschilling &: Marzari 



1996 



Rasio fc Ford 



1996 



Ford et al.lll999l ) suggested that three-body interaction between a cen- 



tral star and two or more planets can lead to the ejection of one of the planets, which takes away 
energy from the system and leaves the third object on a smaller and eccentric orbit, which is finally 
circularized through tidal interaction between the planet and its host star. iMurray et al.l (jl998l ) 
suggested that resonance interactions between a planet and a disk of planetesimals can also lead 
to inward orbital migration of the planet. Another plausible explanation is provided by the theory 
of the gravitational interaction between a protoplanetary disk and an embedded planet, which was 
formulated, as mentioned in the beginning of this article, more than a decade before the discovery 



^see http : / /exoplanet .eu / , 
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of extrasolar planets. 



In the third mechanism, a disk interacts with an embedded planet through their mutual grav- 
itational forces. The planet causes the formation of spiral waves inside the disk at the Lindblad 
resonances; as a result of the density asymmetry induced by spiral waves, the inner disk exerts 
a positive gravitational torque onto the planet and the outer disk exerts a negative gravitational 
torque onto the planet. According to analytical analy s is, th e overall torque is generally negative 
and, hence, forces the planet to migrate inward. IWardI (|l997l ) pointed out that two major kinds of 
migration exist based on the mass of the embedded protoplanet. In the so called Type I migration, 
the planet's mass is small and the response of the disk is linear; the migration speed is very fast so 
that the migration timescale is as short as ~ 10^ yr for a planetary core of lOM® in a sufficiently 
viscous disk. In Type II migration, the protoplanet has enough mass to open a gap inside the disk, 
and migrates on a much longer viscous timescale. Besides classical analytical analysis, many groups 
have studied t he nonlinear evolution of a disk-planet system using numerical multi - dimensional hy- 
drodvnamics (|K 



2003 



Bate et al 



my 



19991; 



Nelson et al. 



2000 



Papaloizou. Nelson. &: Masset 



20031 . and references there in). Recently, 



2001 



de Val-Borro et al 



.Nelson fc Willvl 
"(j2006l ) also carried 

out a collaborative multi-code comparison research. These numerical simulations showed that the 
nonlinear evolution of the orbital migration of a planet inside a disk agrees with linear analysis in 
a qualitative manner. 

However, nature is always not as simple as human bein gs have naive ly th ought. Unexpected 
pheno mena often arose in the non-linear regime. For example, iMassetl ()200ll ) and lMasset &: Papaloizou 
(| 20031 ) suggested that the corotation torque originated from the coorbital region may play a very 
different role from that of Lindblad torques; this leads to a third kind of migration often referred 
to as Type III migration or runaway migration, in which the migration ha ppens on a timescale as 



short as a few tens of orbits and can be directed outward in some cases. Klahr fc Bodenheimer 



(|2003l ) described a Baroclinic Instability in non-barotropic disks that may contribute to vorticity 
and global trubulences; Th ey argued that strong vortici ti es may contribute t o the rapid for r nation 
of Jupiter-size gas planet (IKlahr &: Bodenheimen l2006l ) . iKoller et al.l ([2003) and iLi et al.l (120051 ) 
showed that the so-called "Rossby-wave instabilities" may develop at the local minima of potential 
vorticity (PV), or vortensity (defined as the ratio between local vorticity and surface density), in 
an inviscid disk with initially uniform vortensity distribution. 



In simulations of lLi et al.l (|2005l ). these non-axisymmetric instabilities lead to the formation of 
vorticies and density blobs, which exert stronger torque onto the planet when they travel around 
its vicinity and bring large oscillations to the total torque acted on the planet. They further argued 
that this mechanism may be possible to change the direction of the migratio n. Non-axisymmetric 
Rossby-wave instabilities are also relevant to the evolution of a si ngle disk (ILovelace et al.lll999l : 
Li et al.ll2000l ) and stellar models with strong differential rotation (lOu &: Tohlind 120061 ). The time 
scale for RWIs to fully develop is generally on tens of dynamical time of a disk (~ 100 orbits), hence, 
they can not be ignored in inviscid or low viscosity cases. Therefore, it worths more labor to study 
their comprehensive role on planetary orbital migration and the robustness of this mechanism. 
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In this paper, we present results from simulations of disk-planet interaction. Our purpose is 
two-folded: first w e wish to provide further comparisons to results presented in the collaborative 
comparison work of lde Val-Borro et alj (|2006l ): on the other hand, we focus more on the development 
of non-axisymmetric RWIs in disk models with non-barotropic equation of state (E OS) and non- 



unifor m in itial vortensity distribution, as in contrast to disk models presented in iKoller et al 



(|2003l ) and lLi et alJ (120051 ) . which are isothermal disks with uniform initial vortensity everywhere. 
In particular, we find out that the PV generation is more effective an d commo n in di sk models with 
non-barotropic EOS. Our results, together with those presented by iLi et al.l (j2005l ). suggest that 
the development of RWIs is very robust in disk-planet systems with different mass distribution, 
rotational laws and EOSs. We also address the impact of RWIs on planet orbital migration and 
formation of close-in super Earths near a giant planet. In section 2, basic equations, methods and 
initial models are described; results from test runs are given in section 3; section 4 concentrates on 
the development of RWIs under certain circumstances; we close with conclusions and discussions 
in section 5. 



2. Basic Equations, Methods and Initial Setup 



To study the interaction between a disk and an emb edded planet requir es coupling hydrody- 
namics and orbital dynamics together. Here, we follow iNelson et al.l (120001 ) and many previous 
investigations to reduce the problem to a two-dimensional (2D) one since the disk is considered 
to be very thin. Three dimensional (3D) investigations will be postponed to future. The fluid 
motion inside the disk is described by the vertically integrated continuity equation ([!]), radial and 
azimuthal components of the Navier-Stokes equation ([2|) and 



as 

dt 



+ V • i^v) 



dt 
dt 



+ V • (J^v^v) 







dP j^^^ ^ f 
dr dr 



IdP 

r d(h 



(1) 

(2) 
(3) 



where S is disk surface density, v is two fluid velocities, P is vertically integrated pressure, fr and 
/(^ are two components of viscous forces, and $ i s the gravitat i onal p otential felt by fluid elements. 
Details regarding viscous terms can be found in iNelson et al.l (j2000l ). The EOS of the disk fluid is 
considered as locally isothermal (INelson et al.ll2000l ) as given in below 



P 



(4) 



where the local isothermal sound speed is = ^^JGM^,/r with disk aspect ratio H/r = 0.05. (As 
will be discussed in details later, this radial variation of sound speed generates vorticity wherever 
an azimuthal density gradient exists.) 
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We further simplify our study to non-self- gravitational systems, in which the self- gravity of the 
fluid is not taken into account for the fluid motion; hence, $ = <&=K-|-$p, where is the potential field 
of the central star and $p is the potential field of the planet, which is given by $p = —Mp/y/r"^ + e-^, 
where Mp is the planet mass and e is taken to be 0.2 of the Roche Lobe of the planet. The initial 
disk model has Keplerian rotational profile and uniform density, which results in an initial radial 

3 

vortensity profile ^(r) that is proportional to r~2 (vortensity is defined as the ratio between local 
vo rticity and surface density ) . The value of density and viscosity are chosen to follow those specified 



m 



de Val-Borro et al.l (120061 ). The disk and the non-rotating coordinate system are centered at the 



central star instead of the center of mass of the system, which locates slightly off the central star. 
To handle the hydrodynamics part, we adopted a legacy code develope d by the astrophysica l grou p 



at Louisiana State University to study star formation ( Tohline 1980 : Williams &: Tohline 198^) 



mass transferring binary stars (jMotl. Tohline. &: Frankll2002l ). and rotating instabilities in neutron 
stars ljOu k. Tohlind 12000 ) . The code is explicit and 2nd order in both space and time. It splits 
the source term and advection term in a manner similar to Zeus (jStone &: Normanlll992l ). Other 
features implemented include Van Leer upwind scheme, artificial viscosity to handle shock, and, 
staggered cylindrical grids. The code is originally three-dimensional, but adapted to 2D in this work. 
At the boundary of our computational grids, mass is allowed to fiow off the gr ids but no infiow is 



allowe d. We also implemented the wave-killing boundary condition specified in lde Val-Borro et al 
(j2006l ) for comparison purpose. 



The equation of motion for the planet is: 



(frp 



G{M, + Mp) 



G 



S(r') 



-{fp — r')r'dr't 



(5) 



{\fp — r'p -|- e^) 2 

where the first term on the right is derived from the relative motion of the planet to the central star, 
the second term is the gravitational force exerted on the planet by the disk. For comparisons with 



previous investigations (|Li et al. 



2005 



de Val-Borro et al.ll2006l ). we turned off the disk's attraction 



in many of our simulations so that the planet stays on a fixed circular orbit at 1 AU ; for certain 
runs, the disk potential was turned on to allow the planet to move so that we can assess RWI's role 
on the migration. A 4th order Runge-Kutta integrator is used to integrate forward the equation of 
motion of the planet during each Courant timestep as requi red by the hydro pa r t. To compute the 
torque acted on the planet, we followed the specification in Ide Val-Borro et al.l (|2006l ). 



(r' 



12 _^ g2)2 



r'dr'd(j) , 



(6) 



(k'-'Pi 

torques from inside and outside the Roche-lobe of the planet are computed separately for both 
inner disk and outer disk. In the following sections, we show time evolution of the total torque 
including contributions from within and outside the Roche-lobe, the temporal behavior of the total 
torque does not change qualitatively when torques from within the Roche-lobe are excluded. 

In equations of motion for both the planet and disk, we neglected the indirect terms due to 
the acceleration to our coordinate system by the disk (the planet's acceleration to the coordinate 
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system is neglected for the disk evolution as well). Since the difference between the center of the 
mass and our coordinate origin is fairly small, we expect the effect of indirect terms won't be 
significant to change our results qualitatively; hence we dropped these terms from our equations 
to save computing time. (The evaluation of these indirect terms involves integrals over the mass 
distribution and is therefore relatively expensive.) Finally, the accretion onto the planet is neglected 
as well. 

The units adopted are the following: the gravitational constant G = 1, length unit is 1 AU, 
and Mi, + M p = 1. These setups help u s to have a direct comparison between our results and those 
presented in lde Val-Borro et al.l ( 20061 ). A Jupiter mass is defined as 10"'^ and a Neptun e mass is 
defined as 10~^. The planet's mass is turned on gradually in the first 5 orbits. As shown in lLi et al 



(j2005l ). very high resolution is required to resolve RWIs. This limits our ability to push simulations 
to much longer time scale. In general, we only advance our simulations to ~ 200 orbits, which is 
sufficient for RWIs to fully develop. 



3. Test Runs for Jupiter mass planet 



In this section, we present results for interaction between a disk and a Jupiter mass planet on 
a fixed orbit as a calibration to our code. According to previous analytical analys is and numerical 



nonlinear studies (jLin Sz Papaloizou 



1986 



Trilling et al. 



1998 



Nelson et al.ll2000l ). a Jupiter mass 



planet is capable of clearing most of the materials around its orbit and opening a fairly wide and 
deep gap in a viscous/inviscid disk. In order to test whether our code can simulate this correctly, 
we have run this configuration for both viscous and inviscid disks on grids with two different 
resolutions, 128 x 384 and 256 x 512, to check for convergence. To avoid lengthy descriptions, we 
only show results for a viscous disk. 

Figure [1] displays radial profiles of surface density S and vortensity ^ averaged over the az- 
imuthal direction at t ~ 0, 100, 200 orbits for the higher resolution run. At the final age, a wide 
and deep gap has been opened from the initial flat density profile, the contrast between the lowest 
and highest density spans over around two orders of magnitude. Inside the gap where the surface 
density is extremely low, the vortensity has increased to a very high level. (Notice that the vorten- 
sity is normalized to its initial value at r = 1.) Two vortensity minima formed around the inner 
and outer edge of the gap, where RWIs are expected to grow. But since surface density inside the 
gap is extremely low, we don't expect it to have significant effect on the evolution. (We note that 
th e existence of viscosity als o affects the development of RWIs.) To compare with plots presented 
in Ide Val-Borro et al.l (|2006l ) , a logarithmic color map of surface density distribution of the disk at 
t ~ 200 orbits is shown in Figure [2j It is observed that spiral arms originating from the planet 
forms in both the inner and outer disks. These two spiral arms are well-known Lindblad waves 
caused by the potential of th e planet and, in turn, exert a negative total torque o nto the planet as 
predicted by classical theory (iGoldreich fc Tremaindll978l : iLin Sz PaDaloizoulll980l). T he structures 
shown in this plot agree qualitatively with those shown in Ide Val-Borro et al.l (|2006l . see their Fig. 
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10). 



Figure [3] illustrates the evolution of the total torque (per unit mass) exerted on the planet by 
the disk for both lower and high er resolution runs. The data have been smoothed using a moving 
box method over 10 orbit period (iNelson et al.ll2000l ). Their behaviors agree closely with each other, 
which suggests that we have achieved convergence on these results. As expected, the total torque 
has a slightly negative value and will drive the planet to migrate inward gradually. The magnitude 
of the torque falls in the s ame range of those presented in a recent collaborative comparison study 
(jde Val-Borro et al.l l2006l . see their Table 6). We also monitored the mass in the computational 
grids, there are around 93% of the initial mass left at the end of the higher resolution run. 

In general, our test runs for Jupiter mass planet agree with previous investigations in a quali- 
tative manner. In the following, we will switch to the situation of a Neptune-mass planet embedded 
in a protoplanetary disk and study the development of RWIs in such configurations. 



4. High Resolution runs for RWI in protoplanetary disk with a Neptune-mass 

planet 

4.1. Results for simulations with a Neptune mass planet on a fixed orbit 

In this subsection, we present results from a series of simulations in wh i ch a Neptune-mass 



plane t is on a fixed circular orbit. Previous investigations (jNelson et al.l l2000l : Ide Val-Borro et al 



20061 ) have shown that a Neptune- ma ss planet can only creat es a sha l lower and narrower gap 
inside a disk due to its smaller mass. iKoller et al.l (120031 ) and iLi et al.l (120051 ) further suggested 
that, in an inviscid disk with a Neptunian mass planet, radial vortensity minima are important 
because they may trigger RWIs to develop. As a result of RWI, non-axisymmetric density blobs 
form and orbit around the central star; the torque exerted on the planet by the disk can undergo 
large amplitude oscillations as these density blo bs pass by the vicinity of the planet and apply 
gravitational perturbations to it. As suggested bv lLi et al.l (|2005l ) . this mechanism may potentially 
slow down the inward migration of a planet and, in some extremely cases, even change the sign of 
the total torque exerted on the planet and, thus, cause it to migrate outward. 



Since the disk models studied in iKoller et al.l ([2003) and iLi et al.l (|2005l ) have an isothermal 
EOS and an initial flat vortensity profile, we wish to investigate if the same mechanism will also 
work in disks with non-barotropic EOS an d non-uniform vorte n sity. With this in mind, we adopted 
Keplerian disk models similar to those of Ide Val-Borro et al.l (j2006l ). which has a locally isother- 
mal EOS and an initial uniform surface density and, hence, an initial monotonically decreasing 
vortensity profile as proceeding outward radially. We carried out simulations of this configuration 
on grids with resolutions of 128 x 384, 400 x 800, and 400 x 1600; a simulation with a resolution 
of 800 X 3200 was also performed for ~ 70 orbits to check for the convergence of vortensity profile. 
The viscosity is turned off in these simulations. 
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Fig. 1. — The top and bottom panels show, respectively, and ^(R) of a viscous disk with a 

Jupiter mass planet averaged azimuthally at t ~ (solid line), 100 (dotted line), and 200 (dashed 
line) orbits for the higher resolution run (256 x 512). 
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Fig. 2. — Logarithmic color map of surface density distribution for a viscous disk with a Jupiter 
mass planet at t ~ 200 orbits for the higher resolution run (256 x 512). The horizontal axis is the 
radial axis ranging from 0.4 to 2.5, and the vertical axis is the azimuthal direction spanning over 
27r range with the planet shifted to around the middle. A letter "P" is labeled next to the location 
of the planet. The color bar represents relative rather than absolute values. 
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Fig. 3. — Time evolutions of the total torque (per unit mass) exerted on a Jupiter mass planet by 
the disk for lower (solid line) and higher (dotted line) resolutions 
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In general, our simulations show that the evolution of the system consists of two stages: the first 
stage is the formation of two spiral arms which exert a negative torque to the planet; the second 
stage is the development of RWIs in valleys (minima) of the radial vortensity profile some time 
after the formation of spiral arms. However, fine structures are only resolved accurately i n higher 



resolution runs, mainly because interactions happen mostly around the co-orbital region (jLi et al 



20051 ). In the following, we show detailed results from a run with a resolution of 400 x 1600. 



Figure H] displays radial profiles of S and ^ averaged over the azimuthal direction at t ~ 
0, 60, 120 orbits. In order to give a better view of fine structures around the co-orbital region, only 
regions from r = 0.6 to r = 1.4 are shown in the plot. Compared to the case of a Jupiter mass 
planet, a Neptune mass planet only opens a relatively shallower and narrower density gap. The 
contrast between the lowest and highest density is around a factor of two. Since the surface density 
inside the gap is not quite low, the vortensity there does not increase to much high level. (Notice 
that the vortensity is normalized to its initial value at r = 1.) However, three vortensity valleys 
formed across this radial region, with one centered around the planet's orbit. Because the density 
inside the gap is still significant, RWI s resulted from these vortensity minima will have important 
effect on the evolution of the system (jLi et al.ll2005l ). 



A linear color map of surface density distribution on a polar plot (iLi et al.ll2005l ) at t ~ 150 
orbits is illustrated in Figure [H Focusing on the formation of high density areas (red/brown 
regions), we observe not only red Lindblad spiral arms, but also, other non-axisymmetric high 
density structures at different locations: the inner edge of the outer disk (brown arc-shaped region 
centered around 7 oclock) , the outer edge of the inner disk (red/brown arc-shaped region around its 
edge); interestingly, the density inside the gap is no longer axisymmetric any more, as suggested by 
the light blue region around 10 oclock. The radial locations of these non-axisymmetric structures 
match exactly with the local mini ma of the vort ensity distribution (see bottom portion of Figure 
HI) as predicted by linear analysis (ILi et al.ll200d ). 



Figure [6] shows the time evolution of the total torque on the planet, the torque from the inner 
and outer disk. The behavior of the total torque is totally different from that in the Jupiter mass 
case: very larg e amplitude oscillations present through most of the simulation. As pointed out by 



Li et al.l (|2005l ). these oscillations are caused by the gravitational pull from those orbiting density 



blobs; each time when they pass the vicinity of the planet, they exert a stronger torque onto the 
planet, the sign of their torques depends on their relative position to the planet. The existence of 
these oscillations complicates the classical picture of the gravitational torque exerted on the planet 
by the disk. 

Figure El exhibits the radial profile of the specific angular momentum J inside the disk around 
the co-orbital region at different time. The solid line denotes for initial Keplerian profile, dotted, 
dashed, and dotted-dashed lines denote for profiles at later times. As the evolution goes on, the 
inner disk exerts a positive torque to the planet and, hence, is giving away angular momentum 
to the planet; whereas, the outer disk exerts a negative angular momentum to the disk, thus, is 
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gaining angular momentum from the planet; therefore, the angular momentum is being transferred 
from the inner disk to the outer disk due to the interaction between the planet and the disk. This 
mechanism results in a dip in the inner disk's J profile but a bump in the outer disk's J profile. 
As a result, an S-shaped structure formed in the vicinity of the planet's fixed orbit. Again, for 
comparison to previous studies, there are around 98% masses left in the grid at the end of the 
highest resolution run. 

In summary, our results of the development of R WIs in an invis cid disk with a Neptune-mass 
planet on a fixed orbit are in agreement with those of lLi et al.l (120051 ) ; how ever, they differ in some 
key physical quantities. In particular, the PV maxima in iLi et al.l (|2005l . see their Figure 2 and 
9) are located at radii outside the coorbital region where shocks present; whereas, PV maxima in 
our simulations reside within the coorbital region around 0.95 < r < 1.05, where shocks are not 
expected to occur. This discrepancy raises some concerns. Because based on the evolution equation 
of vorticity C,: 

dC VS X VP 

the vorticity can only be generated at regions where pressure gradient and density gradient do not 
align with each othe r; such generati on of vorticity normally happens around shocks in barotropic 
disks as illustrated in lLi et ah (2005), which contributes to vortensity maxima outside the coorbital 
region. 

To understand the mechanism that generates vortensity (PV) maxima within the coorbital 
region in our simulations, we note that our disk models have a non-barotropic EOS with a radial 
variation of sound speed Cg, which also brings in misalignment of pressure gradient and density 
gradient wherever azimuthal density gradient appears (Vc^ x VS ^ 0). This misalignment acts 
as a source term for the generation of vorticity within the coorbital region, where no shock but 
strong azimuthal density gradient presents. Therefore, the large density depre ssion in the same 
region naturally gives birth to vortensity increment. According to lLi et al.l (|2005l ). lower resolution 
may result in unphysical generation of vortensity. To further test that if our results are resolution- 
dependent, we carried out a run with resolution 800 x 3200 for 70 orbits (due to the limitation 
of computational resources). Figure [8] shows the vortensity profile around coorbital regions for 
runs with different resolutions at t ~ 70 orbits. They all exhibit vortensity maxima within the 
coorbital region, which indicates that our results are resolution- independent. As a further evidence 
that vortensity is generated around where azimuthal density gradient exists. Figure [9] and [10] il- 
lustrate the distribution of vortensity increment and azimuthal density gradient, respectively, for 
the run with resolution 800 x 3200 at t ~ 70 orbits. For a better view, they are zoomed in the 
neighborhood of the planet. It is observed that vortensity is generated within the Roche lobe, 
where strong azimuthal density gradient, instead of shocks, exists. Such a generating mechanism 
of vortensity inside a planet's R oche lobe has the same origin as baroclinic instability discussed 



m 



Klahr &: Bodenheimed (|2003l ). which contributes to global turbulence within the disk. On the 
other hand, we also note that the oscillations of the torque acted on a neptu ne mass planet appear 



much earlier in our simulations (~ 50 orbits in our Figure 6) than they do in lLi et al 



([20051, ~ 200 
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orbits in their Figure 6), which suggests that PV (vortensity) generation is even more effective and 
common in non-barotropic disk models. 

As show in Figure [TOl the azimuthal density gradient flips the sign across the planet; thus, we 
expect the vortensity increment also flip the sign there. However, the vortensity change within the 
Roche-lobe shown in Figure [9] is always positive. We don't fully understand the underlying physics 
of this phenomenon, but suspect that this is possiblely due to vortensity mixture along librating 
stream lines in the horse region, i.e., where fluid elements make the U-turn. A detailed and vigrous 
analysis of fluid motion and vortensity mixture around a Neptune-mass planet embedded in a non- 
barotropic disk is beyond the scope of this paper, and will be studied in a separate investigation. 

Since the generation of vortensity in non-barotropic disks turns out to be quite universal, the 
second instability (RWIs) triggered around vortensity extrema is probably more common in realistic 
situations than what has been thought previously; it will have strong influence on the migration 
of a neptune mass planet in an inviscid disk. Next, we will discuss the case in which the planet is 
allowed to move freely under the torque action from the disk. 

4.2. Simulation for a freely moving Neptune-mass planet 

In order to assess possible roles of RWIs on the migration of a planet, we carried out two simu- 
lations in which a Neptune mass planet is allowed to move under the influence of the gravitational 
torque from the disk. The resolutions are 400 x 800 and 400 x 1600. In these two simulations, the 
planet is fixed on the initial circular orbit for the first 50 orbits (slightly longer than the libration 
period for the neptune mass planet) to allow the disk to respond properly, it is then released to 
move freely under the influence of the torque from the disk. The results from these two simulations 
are similar to each other, which suggests that convergence is achieved. We use the higher resolution 
run to present the evolution of the system. 

Figure [TT] illustrates radial profiles of S and averaged over the azimuthal direction at t ~ 
50, 100, 150, 200 orbits. It is observed that at t = 50 orbits the planet has created a shallow gap 
within the disk around r = 1, the density level is enhanced at both the inner and outer edge of the 
gap (r ~ 0.88 and r ~ 1.12). Then, the planet migrates inward when t > 50 orbits as an effect 
of an overall negative total torque. As it moves in, the surface density gap and vortensity valleys 
also migrate with it. On the other hand, a shallower density gap and a vortensity valley are still 
preserved around its initial orbital location. 

Figure [12] exhibits the time evolution of the total torque acted on the planet and its radius 
to the central star. The overall torque exerted on the planet is negative in average throughout 
the simulation, but undergoes some changes. After the whole system settles down in the first 50 
orbits, large amplitude oscillations with certain period show up; the torque becomes increasingly 
negative when the planet starts to migrate inward; it then reaches its peak in magnitude at t ~ 100 
orbits (~ 50 orbits after releasing); after that, it is weakened and then returns back to a very small 
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Fig. 4. The top and bottom panels show, respectively, and ^(i?) of the disk with a Neptune 

mass planet averaged azimuthally at i « (solid line), 60 (dotted line), and 120 (dashed line) 
orbits. 
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Fig. 5. — The same as Fig. [2] but in linear scale for an inviscid disk with a Neptune mass planet at 
t ~ 150 orbits. A letter "P" is labeled next to the location of the planet. The color bar represents 
relative rather than absolute values. 
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Fig. 6. — Time evolutions of total torque (top), torques from the inner (middle) and outer (bottom) 
disks exerted on a Neptune mass planet for the highest resolution run (400 x 1600). Raw data is 
shown. 
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Fig. 7. — Radial profiles of specific angular momentum around the co-orbital region of the disk with 
a Neptune mass planet at different times. The solid line denotes for initial Keplerian profile, dotted, 
dashed, and dottcd-dashed lines denote for profiles at later times. Data have been normalized to 
the initial value at r = 1. As evolutions goes on, angular momentum is being removed in the inner 
region but added to the exterior region due to the interaction between the disk and the planet. 
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Fig. 8. — Vortensity profiles around the coorbital regions for runs with different resolutions at 
t ~ 70 orbits. All of them show maxima within the coorbital region. 



Fig. 9. — Distribution of vortensity increment around the planet at t ~ 70 orbits for the run with 
resolution 800 x 3200. (The initial backrougnd vortensity is subtracted.) Shocks are well resolved. 
Vortensity is generated within the coorbital region, especially in the Roche lobe, where no shock 
exists. The azimuthal density gradient within the same region acts as a source term of vortencity 
generation (see Figure [TO]) . Note that the vortensity increment within the Roche- lobe is always 
positive, which contradicts our expectation. We suespect this may be resulted from vortensity 
mixture in this region. 




Fig. 10. — Distribution of azimuthal density gradient at the same time as that of Figure [9l Strong 
density gradient is observed within the Roche lobe and coorbital region. 
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negative value with smaller oscillations that have a different period from that of earlier oscillations. 

We carefully analyzed the evolution between t = 80 ~ 120 orbits and found that the sudden 
change of the torque's magnitude is due to the development of a relatively large density blob 
between the planet and the edge of inner disk as a result of RWIs. During this time, the density 
blob accumulates materials and augments very quickly because the planet is penetrating the high 
density region near r ~ 0.88 (which is introduced by the planet's potential during the first 50 orbits, 
see the density profile at i = 50 orbits). Consequently, the density blob exerts a stronger torque 
onto the planet causing the oscillation amplitude of the total torque to increase, the averaged torque 
also becomes increasingly negative during t = 80 ~ 100 orbits. After the planet passed the high 
density region at f ~ 100 orbits, the large density blob saturates and wanes; the torque gradually 
returns back to a value close to zero. However, smaller density blobs do survive and contribute to 
smaller oscillations of the total torque in the late evolution. Figure [13] shows the density colormaps 
of the disk at t ~ 100 orbits when the planet is penetrating the high density region, a pronounced 
density blob (green arc-like structure) in the edge of the inner disk is clearly visible. For comparison, 
similar plot for t ~ 200 orbits is shown in Figure [T^ where previous large density blob diminishes 
and only small density blob survives both in the inner disk and outer disk. In order to have a 
better view of the time behavior of the density azimuthal asymmetry between the planet and the 
outer edge of the inner disk, the radially-averaged density azimuthal distributions of this ring-belt 
at different times are ploted in Figure [151 Besides some sharp spikes that signal the locations of 
the planet, density azimuthal asymmetry with different degrees is clearly seen. We note that the 
density blob has the highest value at t ~ 104 orbits, which coincides with the time when the torque 
has the negative maximum value and when the planet is penetrating the high density region at 
r ~ 0.88. This further supports our analysis on the correlation between the torque magnitude and 
the augmentation of the density blob. 

As a consequence of this interesting phenomenon, the planet's radius decays quickly during 
t = 80 ~ 120 orbits, but the decaying slows down after t « 110 orbits (see bottom panel of Figure 
[12] for this non-monotonic behavior). We measured the migration timescale Tmig, which is defined 
as the time needed for the planet to migrate from r = 1 to the central star given a drifting speed 
measured at the end of our simulation. Our measurement yields Tmig ~ 1300 orb its, which is almost 



an order of magnitude lower than that of type II migration from a similar study (iNelson et al.ll2000l . 
see their Table 1. Note that their planet also locates at r = 1 initially.). This timescale is consistent 
with type I migration timescale, which ought to be one or two orders of magnitude lower than that 



of type II (|Wardlll997l l 



We note that there is a trend for the averaged torque to approach zero at late times, this is 
probably because the mass of the disk interacting with the planet decreases as it moves in. On the 
other hand, the oscillation amplitude of the torque induced by RWIs remains relatively constant, it 
is possible that as the planet moves closer to the host star a turning point may eventually appear, 
where the RWI-induced torque beats over the Lindblad torque. If this is true, the planet may be 
halted at certain radial range. Since the planet is very close to the inner edge of our computational 
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grid at the end of this simulation, we decided to stop the simulation and will address these issue in 
future investigations. 

In the limited time evolutions of a freely moving planet presented here, RWIs do not change 
the overall picture of inward migration; but they have significant influence on the torque exerted 
on the planet and consequently make the migration speed of the planet non-monotonic. Further 
investigations are needed to study the influence of other factor that may potentially bring in sig- 
nificant changes, such as the "releasing time" of the planet, self gravity of the disk fluid, etc. One 
important issue is how to properly evaluate the torque from materials near the planet, which is 
embedd ed in the disk . In this paper, we cornputed this bru teforcely in 2D. However, many previous 
studies (jMassetl I2OO2I : iD'Angelo et al.ll2005l : iLi et al.ll2005l ) used quasi-3D treatment, which tends 
to reduce the torque. It is very likely that if proper 3D treatment is adopted, the total torque 
will be reduced and the influence of RWIs will be much stronger. This could potentially cause the 
migration to be reversed. 



Conclusion and Discussion 



We have carried out high resolution simulations on the interaction between a protoplanetary 
disk and an embedded planet with emphasis on detailed interplay between a disk and a planet 
under the influence of baroclinic generation of vortensity and non-axisymmetric RWIs. Our results 
are consistent with classical analysis on the interaction between a protopla netary disk and an 
emb edded pl a net th rough Lindblad torques. We confirmed previous studies bv lKoller et al.l (j2003l ) 
and iLi et al.l (|2005l ) that non-axisymmetric RWI is likely to develop under certain circumstances 
and have an important influence on the migration of a planet inside an inviscid disk. We also 
found that the generation of vortensity (PV) is more common and effective in disks with non- 
barotropic EOS through the baroclinic instability, which further favors the development of RWIs. 
As the asymmetry of the density distribution induced by RWIs becomes prominent, the resulting 
density blobs exert periodical and enhanced gravitational pull onto the planet as they pass by the 
vicinity of the planet, which causes the total torque received by the planet undergo large amplitude 
oscillations. Although in our current simulations, RWIs did not change the overall picture of an 
inward migration, they definitely have very important and interesting effect on the migration speed. 
As a side effect of an inwardly migrating planet, RWIs introduce nonaxisymmetric density blobs 
along the way. These enhanced density blobs with s trong vortices may help rap i d formation of 



new planet cores within them in a way suggested in (iKlahr &: Bodenheimer 



2006 



Petersen et al 



20061 ). especially in inner regions of c ircumstellar disks wher e rapid precipitation and coagulation of 
solid materials are likely to happen ( Silyerstone et al. 20061) . If these new-born cores could survive 
during the migration of a giant planet ( ChambersI 2006|: Ravmond et al. 2006 ). they may produce 
Earth-l ike planets in the Hab itable Zones (jji et al.ll2006l ) or Hot Earths interior to a close-in giant 
planet (jRavmond et al.ll2006l ). 



The study presented here is just a very limited step among many efforts toward understanding 
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Fig. 11. — The same as Fig. |3]but for a run at t 50 (solid line), 100 (dotted line), 150 (dashed 
line), and 200 (dash-dotted line) orbits, in which a Neptune mass planet is allowed to move after 
50 orbits. 
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Fig. 12. — Time evolutions of the total torque (top) exerted on a Neptune mass planet by the disk 
and the radius of the planet (bottom) for the run in which the planet is allowed to move. Raw data 
is shown. A non-monotonic behavior for the migration speed is observed. 
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Fig. 13. — The same as Fig. [5] but at t = 103 orbits for a run, in which a Neptune mass planet is 
allowed to move. The planet is penetrating the high density region around r ~ 0.88 and a fairly 
large density blob forms (green arc-like structure at the edge of inner disk). 
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Fig. 14. — The same as Fig. [13] but at t = 202 orbits. The planet has passed the original high 
density region around r ~ 0.88. Previous large density blob wanes, smaller density blobs survive 
at several locations. 
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Fig. 15. — The radially-averaged azimuthal density distribution of the ring belt between the planet 

and the outer edge of the inner disk from a run in which a ncptune mass planet is allowed to 
move. Note that the density blob has the highest value at ~ 3 and t ~ 104 when the planet 
is penetrating the high density region around r ~ 0.88. (The spikes reveal the locations of the 
planet.) 
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disk-planet interactions. Compared to realistic situations, there are so many physics missing in 
our much-simplified simulations, the influence of RWIs on planet migration is still not clear in 
cases where all the physics are correctly taken into account. For example, if the self-gravity of 
the disk ignored in this study is included in the dynamics of the fluid, it is likely to enhance the 
RWIs because the gravitational potential field is made more asymmetric by the non-axisymmetric 
density distribution and more masses will be trapped into potential wells of these density blobs. 
Enhanced RWIs will exert a stronger positive torque onto the planet and may greatly reduce the 
migration speed. On the other hand, because RWI develops mostly around the coorbital region 
where corotation torques are originated, it may have important effect on type III migration. These 
issues need to be addressed in future investigations. 
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